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Abstract 

As shown in m, under some structural assumptions, working on con¬ 
gested traffic problems in general and increasingly dense networks leads, 
at the limit by T-convergence, to continuous minimization problems posed 
on measures on generalized curves. Here, we show the equivalence with 
another problem that is the variational formulation of an anisotropic, 
degenerate and elliptic PDE. For particular cases, we prove a Sobolev 
regularity result for the minimizers of the minimization problem despite 
the strong degeneracy and anisotropy of the Euler-Lagrange equation of 
the dual. We extend the analysis of [6] to the general case. Finally, we 
use the method presented in [5] to make numerical simulations. 

Keywords: traffic congestion, Wardrop equilibrium, generalized curves, 
anisotropic and degenerate PDEs, augmented Lagrangian. 


1 Introduction 

Researchers in the field of modeling traffic have developed the concept of 
congestion in networks since the early 50’s and the introduction of the 
notion of Wardrop equilibrium (see [22]). Its important popularity is due 
to some applications to road traffic and communication networks. We will 
describe the general congested network model built in ns in the following 
subsection. 

1.1 Presentation of the general discrete model 

Given d e N, d > 2 and H a bounded domain of R d with a Lipschitz 
boundary and e > 0, we take a sequence of finite oriented networks = 
(N e ,E s ) whose characteristic length is e, where N e is the set of nodes 
in and E £ the set of pairs (a:, e) with x £ N £ and e £ I d such that 
the segment [a:, x + e] is included in fi. We will simply identify arcs to 
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pairs ( x,e ). We assume \E E \ = max{|e|, there exists x such that (x,e) £ 
E E } = e. . 

Masses and congestion: Let us denote the traffic flow on the arc 
( x,e ) by m E (x,e). There is a function g E : E e x R+ —» R+ such that for 
each (x,e) £ E E and m > 0, g E (x,e,m ) represents the traveling time of 
arc ( x, e ) when the mass on (x, e) is m. The function g E is positive and 
increasing in its last variable. This describes the congestion effect. We 
will denote the collection of all arc-masses m e (x, e) by m £ . 

Marginals: There is a distribution of sources ft = Y2 xeN e /- {x)5 x 
and sinks /+ = Y2 xeN e f+(x)5 x which are discrete measures with same 
total mass on the set of nodes N E (that we can assume to be 1 as a 
normalization) 

£ /-(*) = £ f+(v) = i- 

x£N e y£N e 

The numbers ft(x) and f+(x) are nonnegative for every x £ N E . 

Paths and equilibria: A path is a finite set of successive arcs (x, e) £ 

E E on the network. C e is the finite set of loop-free paths on and may 
be partitioned as 

= U = U = U ( 

(x,y)eN c xN e xeN e 

where . (respectively C £ tV ) is the set of loop-free paths starting at the 
origin x (respectively stopping at the terminal point y ) and Ct, y is the 
intersection of CJ r and Cf y . Then the travel time of a path 7 £ C e is 
given by: 

Tme( 7 )= £ g e (x,e,m £ (x,e)). 

(x,e)C7 

The mass commuting on the path 7 £ C £ will be denoted w e (y). The 
collection of all path-masses w e ('y) will be denoted w s . We may define 
an equilibrium that satisfies optimality requirements compatible with the 
distribution of sources and sinks and such that all paths used minimize 
the traveling time between their extremities, taking into account the con¬ 
gestion effects. In other words, we have to impose mass conservation 
conditions that relate arc-masses, path-masses and the data ft and /+: 

/-(*):= £ w e ( 7 ), f+(y) := £ w e (-y), V(®, y) £ N e x N e (1) 
7SCI,. 7 eC‘ y 

and 

m E (x,e)= £ w e (-y),V(x,e) £ E £ . (2) 

7€C e :(a;,e)C7 

We define Tfje to be the minimal length functional, that is: 

T £ e(x,y):= min £ g e (x, e, m £ (x, e). 

,y (x,e)C7 

Let U(ftJ+) be the set of discrete transport plans between ft and /+, 
that is, the set of collection of nonnegative elements (<fi £ (x, y))( X:y )eN^ 2 
such that 

£ <p e (%,y) = /-(*) and £ <p e {x,y) = f+(x), for every (x,y) £ N E xN e . 

yeN £ xeN £ 
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This results in the concept of Wardrop equilibrium that is defined precisely 
as follows: 

Definition 1.1. A Wardrop equilibrium is a configuration of nonnegative 
arc-masses nT : (x,e) —> ( m s (x,e )) and of nonnegative path-masses w e : 

7 —> w e ('y), that satisfy the mass conservation conditions 0 and |2]| and 
such that: 

1 . For every ( x,y) £ N s x N e and every 7 £ C% y , ifw e ( 7) > 0 then 

Tmc( 7 )= min r^ e ( 7'), ( 3 ) 

T'eoS,„ 

2 . if we define n E (cc, y) = w e ( 7 ) then IF is a minimizer of 

eJjt f e y) • ( 4 ) 

Condition 0 means that users behave rationally and always use short¬ 
est paths, taking in consideration congestion, that is, travel times increase 
with the flow. In [MS], the main discrete model studied is short-term, 
that is, the transport plan is prescribed. Here we work with a long-term 
variant as in 00 ■ It means that we have fixed only the marginals (that 
are ft and /+). So the transport plan now is an unknown and must be 
determined by some additional optimality condition that is 0 . Condi¬ 
tion Q requires that there is an optimal transport plan between the fixed 
marginals for the transport cost induced by the congested metric. So we 
also have an optimal transportation problem. 

1.2 Assumptions and preliminary results 

A few years after the work of Wardrop, Beckmann, McGuire and Winsten 
[2] observed that Wardrop equilibria coincide with the minimizers of a 
convex optimization problem: 

Theorem 1 . 1 . A flow configuration (w e ,nF) is a Wardrop equilibrium 
if and only if it minimizes 

nm 

G e (x, e, m e (x, e)) where G s (x, e, m) := / g e (x,e,a)da (5) 

(:c,e)6B e "' 0 

subject to nonnegativity constraints and the mass conservation conditions 

0 -0. 

The problem 0 is interesting since it easily implies existence results 
and numerical schemes. However, it requires knowing the whole path flow 
configuration w e so that it may quickly be untractable for dense networks. 
However a similar issue was recently studied in m- Under structural 
assumptions, it is shown that we may pass to a continuous limit which 
will simplify the structure. Here, we will not see all these hypothesis, only 
the main ones. So we refer to na for more details. 

Assumption 1. The discrete measures {e ^ 1 fl) e >o and (ei _1 /I) £ >o 
weakly star converge to some probability measures /_ and f+ on SI ; 

lim £ d/2-1 V {<p(x)ft(x)+if(x)fl{x)) = f <pdf-+f ipdf + y(tp,ip) £ C'(SI) 2 . 

£ ^ 0+ -Five J n Jn 
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Assumption 2. There exists N £ N, {'U; s }fc=i,...,jv £ C 1 (R d ,S d_1 ) iv and 
{ck}k=i,...,N £ C 1 (n,K+) JV such that E E weakly converges in the sense 
that 


lim 

E—>0 + 


E 

( x,e)£E £ 


N'V 


H 


QxS d_1 


(p(x,v)0(dx,dv),\/ip G C(f 2 xS d 


where 0 G + x 1 ) and 6 is of the form 

N 

9(dx, dv) = y^ c k (x)5 Vk{x) dx. 
k= 1 


Moreover, there exists a constant C > 0 such that for every (x,z,£) £ 
R d x § d_1 x R+, there exists Z £ R+ such that \Z\ < C and 


Z ■ £ = min 


Z-^Z 


N 

(zi, ..., ZN ) £ R+ and ^ z k v k (x) 
k =1 



( 6 ) 


The Ck’s are the volume coefficients and the v k ’s are the directions in 
the network. The next assumption focuses on the congestion functions g £ . 

Assumption 3. g s is of the form. 


g E {x,e,m) = |e| d/2 g , Ve > 0, (x,e) £ E E ,m > 0 (7) 

where g : Q. x S d_1 x R+ H > R is a given continuous , nonnegative function 
that is increasing in its last variable. 

We then have 

G e {x, e, m) = \c\ d G lx,-^-r, , ^ /2 ) where G(x, v, m) := / g(x,v,a)da. 

\ \ e \ \ e \ J Jo 

We also add assumptions on G: 

Assumption 4. There exists a closed neighborhood U of Q such that 
for k = 1,..., N, Vk may be extended on U in a function C 1 (still de¬ 
noted Vk). Moreover, each function ( x,m ) £ U x R+ H > G(x,Vk(x),m) is 
Caratheodory, convex nondecreasing in its second argument with G(x , Vk{x), 0) 
0 a.e. x £ U and there exists 1 < q < d/(d — 1) and two constants 
0 < A < A such that for every (x , m) £ U x R+ one has 


A (m q — 1) < G(x,v,m ) < A (m q + 1). ( 8 ) 


The q-growth is natural since we want to work in L q in the continuous 
limit. The condition on q has a technical reason. It means that the con¬ 
jugate exponent p of q is > d, which allows us to use Morrey’s inequality 
in the proof of the convergence m)- The extension on U will serve to 
use regularization by convolution and Moser’s flow argument. Examples 
of models that satisfy these assumptions are regular decompositions. In 
two-dimensional networks, there exists three different regular decomposi¬ 
tions: cartesian, triangular and hexagonal. In these models, the length 
of an arc in E E is e. The CkS and v k ’s are constant. In the cartesian 
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case, N = A, (vi,V 2 ,v 3 ,V 4 ) := ((1,0), (0,1), (-1, 0), (0,-1)) and c k = 1 
for k = 1,..., 4. For more details, see [Ill- 

Now, before presenting the continuous limit problem, let us set some 
notations. 

Let us write the set of generalized curves 

£ = {(7 ,P) ■■ 7 G ^’“([O, l],n),p € V ' 7 n L'QO, 1])"}, 


where 


( N 

V 1 = < p : t e [0,1] -»> p{t) G M+ and 7 (t) = Vjk( 7 (t)) Pfc(t) a.e. 

I k =1 

We can notice that 77, is never empty thanks to [Assumption 2| Let us 
denote Q £ Q(f~, /+) the set of Borel probability measures Q on £ such 
that the mass conservation constraints are satisfied 


Q(f-,f+) := {Q £ M\(C) ■. e 0# Q = = f+} 

where et( 7 , p) = 7 (t),t £ [0,1], ( 7 , p) £ C. For k = 1,. .., N let us then 
define the nonnegative measures on SI x § d_1 , m® by 

[ ip(x,v)dm®(x,v) = [ (( ip^(t),v k (^(t)))p k {t)dt\dQ(^,p), 
JcixS d ~ 1 JC \J0 / 

(9) 

for every <p £ C(fl x S d 1 ,R). Then write simply rrfi = m k ’ non¬ 

negative measure on SI x §“ _1 . Finally assume that 

Q q (f-,f+) := {Q £ Q(f-J+) : m Q £ L q (9)} ? 0. 


It is true when for instance, /+ and /_ are in L q (Q.) and SI is convex. 
Indeed, first for Q £ Af + (W' 1 ,oo ([0,1], SI)), let us define Iq £ _M + (SI) as 
follows 


f <p di Q = [ (( dQ(y) for ip £ C(SI,I 

J n Jw 1 '« > ([o,ii,n) \J 0 / 


It follows from the regularity results of [10, 121] that there exists Q £ 
Ad + (W 1 ’° o ([0,1],SI)) such that eo#Q = /-, = /+ an d iq £ L q . 

For each curve 7 , let p 
the existence due to 


£ 'P'y 


such that ^fePfc(i) — C| 7 (I)I ( we have 
Assumption 2 I. Then we set Q = (id, p )^Q. We 


have Q £ Q q (f-, /+) so that we have proved the existence of such kind of 
measures. 

Then Wardrop equilibria at scale e converge as e —> 0 + to solutions of 
the following problem 


inf 

QeQ q (f~,f+) 


J 

If! 


fixS'*- 1 


G(x , v, m.(x, v))6(dx, dv) 


( 10 ) 


(see na). Nevertheless this problem ( |10| ) is posed over probability mea¬ 
sures on generalized curves and it is not obvious at all that it is simpler to 
solve than the discrete problem ©• So in the present paper, we want to 
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show that problem (101 is equivalent to another problem that will roughly 
amount to solve an elliptic PDE. This problem is 


inf inf 

seP" {J axS d-i 

where 

V° = { g : n x § d 


G(x, v, g(x, v)) 6(dx, dv)-, —div cr = f f , (11) 


Vx G n, a(x) = y^ v k (x)g(x, v k {x)) 


f = f + — /_ and the equation —div(cr) = / is defined by duality: 


f Vm • a = [ udf, for all u G C 1 (f2), 

Jn J n 

so the homogeneous Neumann boundary condition a ■ un = 0 is satisfied 
on dfl in the weak sense. For the sake of clarity, let us define 


g(x,a):= inf c k (x)G(x, v k (x), g k ) ■= inf G(x, g) 

q€.'P% *■—' Q ev% 

k =1 


where 


VZ = \ g e R+; cr = y ^V k (x)g k > and G(x,g) := y^ c k {x)G(x, v k (x), g k ), 


for x G Q, a G 


We recall that the c k ’s are the volume coefficients 


in 8. (J is convex in the second variable (since G is convex in its last 


variable).The minimization problem (111 can then be rewritten as 


inf < [ g(x, u(x)) dx\ — div a = / 1 . (12) 

{Jn J 

This problem ( |12[ ) looks like the ones introduced by Beckmann [3] 
for the design of an efficient commodity transport program. The dual 


problem of (121 takes the form 


sup 

uew 1 -p(n) 


j udf—f G*(x, Vu(x)) dx 
Jn Jn 


(13) 


where p is the conjugate exponent of q and Q* is the Legendre transform 
of g(x,-). In order to solve (12 I, we can first solve the Euler-Lagrange 


equation of its dual formulation and then use the primal-dual optimality 
conditions. Nevertheless, in our typical congestion models, the functions 
G(x,v,-) have a positive derivative at zero (that is g(x,v, 0)). Indeed, 
going at infinite speed - or teleportation - is not possible even when there 


is no congestion. So we have a singularity in the integrand in (121. Then 
G* and the Euler-Lagrange equation of ( | 13| ) are extremely degenerate. 
Moreover, the prototypical equation of [3 is the following 

-div ( ( | VW |-i)r^)=/- 
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Here, for well chosen g , we obtain anisotropic equation of the form 

N d 

-EE di [bk(x)vki(x)(X7u ■ Vk(x) - 8kCk{x)Y_ J. : ] = /. 

fc=i i=i 


where Vk(x) = (vki(x),... ,Vkd(x)) for k = 1,..., N and x £ fi. In the 
cartesian case, we can separate the variables in the sum but in the hexag¬ 
onal one (d = 2), it is impossible. The previous equation degenerates in 
an unbounded set of values of the gradient and its study is delicate, even 
if all the Sk’s are zero. It is more complicated than the one in [6]. Indeed, 
the studied model in [B] is the cartesian one and the prototypical equation 
is 






/• 


The plan of the paper is as follows. In Section 2, we formulate some 
relationship between (101 and (12 I. Section 3 is devoted to optimality 
conditions for 


(jT2J) in terms of solutions of ( |13| ). We also present the 
kind of PDEs that represent realistic anisotropic models of congestion. In 
Section 4, we give some regularity results in the particular case where the 
Ch s and the Vk s are constant. Finally, in Section 5, we describe numerical 
schemes that allow us to approximate the solutions of the PDEs. 


2 Equivalence with Beckmann problem 

Let us study the relationship between problems M and GD- We still 
assume that all spec ified hypothesi s in |Section 1| are satisfied. Let us 
notice that thanks to Assumption 2 for every a £ L q (Q, R d ), there exists 
g £ V 7 such that g £ L q (9) and g minimizes the following problem : 

inf < / G(x,v, g(x,v)) 9(dx,dv) > . 

gSV 7 l,Jn x gd-l J 

For g £ V a , define g : fl —> R+ where Qk(x) = g(x,Vk(x)), for every 
x £ k — 1,..., N. Now, we only consider g that we simply write g (by 
abuse of notations). 

Theorem 2.1. Under all previous assumptions, we have 


inf (jlOb = inf (12). 


Proof. We adapt the proof in [6|. We will show the two inequalities. 
Step 1: inf ( |10| ) > inf (121. 

Let Q £ Q q (f + ,f~). We build £ L 9 (fl,R d ) that will allow us to 
obtain the desired inequality, we define it as follows : 

[ <pda Q =[ [ tpiiit)) ■ y(t)dt dQ(7,p),V^ £ C(f2,R d ). (14) 

Jn JcJ o 

In particular, we have that —div a® = f since Q £ Q(/_,/+). We now 
justify that 

N 

a Q (x)= / vm Q (x,v) dv = V' Vk{x)ufi(x, Vk{x)) a.e. x £ Ul. 

• / s d - 1 Z=i 
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Recall that for every £ G C{Vt x S' 

N 


f £dm® = [ [ ( ^2^('y(t),v k ('y(t)))p k (t) ) dt dQ{^,p). 
JOxS ^- 1 JC Jo \ k=1 J 

By taking £ of the form £(x,v) = tp(x) ■ v with tp £ C(fi,M“), we get 

rl / N 


[ ip(x) ■ vdm®{x,v) = f [ (y2p k (t)ip{'r{t))-v k {'y{t)))dtdQ{y,p) 
J nxS 1 *- 1 Jc Jo \ k=1 J 

= f <p da®. 

Jsi 

Moreover, since m® > 0, we obtain that m® £ P aQ (and so that a® £ L q ) 
and the desired inequality follows. 


Step 2: inf (101 < inf (121. 

Now prove the other inequality. We will use Moser’s flow method (see 
00 ESI) and a classical regularization argument. Fix 5 > 0. Let cr £ 
L q (n,R d ) and g £ V a n L q (n,R N ) such that 


/, 


nxs ^- 1 


G(x , v, g(x, v)) 9(dx, dv) < inf (121 + 8 


with — div a = f. We extend them outside by 0. Let then rj £ C“?°(R d ) 
be a positive function, supported in the unit ball Bi and such that f Rd y = 
1. For e < 1 so that := S2 + eB\ (s U, we define if(x) := e~ d r/(e~ 1 x), 
a e := rf * a and g k {x) := 1 f * g k (x) for k = 1,..., N. By construction, 
we thus have that a e £ C°°(fl E ) and 

— div (<t £ ) = /+ — /! in and cr e = 0 on <9fi e , 

where f± = y E * (/±1q) + e. But the problem is that we do not have 
g e £ V a . We shall build a sequence ( P E ) in V a that converges to p in 
L g (U,R N ). Notice that 

N r 

o- e (x) = ^2 rf{y)Qk{x - y)v k (x - y) dy 

k=1 ^ 

N A T r 

— y^ J QUx)vk{x) + 'y] / y e (y)g k (x-y)(v k (x-y)-v k (x))dy 

k =1 k =1 J 

There exists p% £ L q (Q, e ) such that for every k = 1,..., N, p% > Q,p% —»■ 0 
and for x G H £ , we have 

N p N 

I e {x) = / y s (y)Qk{x - y)(v k {X -y)- v k (x)) dy = ’Y^p e k {x)v k (x). 

k= 1 ' k=l 

Such a family exists since I s £ L q and I e —> 0 (by using t he fact that the 
v k ’s are in C 1 ([/)) and we can estimate p% with I £ due to Assumption 2 


Then if we set P £ = g £ + p £ , we have P £ £ V a and P £ —> g in L q . 






Define g E (t, x) := (1 — t)fl(x) + tf+(x) Vi € [0,1], x £ fl £ , let then X s 
be the flow of the vector field v e := a E /g E , that is, 

f Xf(x)r=V S (t,X!(x)) 

\a'q(x) = x, (t,x) e [0,1] x 

We have dtg E + div (g E v E ) = 0. Since v E is smooth and the initial data 
is gr e (0, ■) = fl, we have Xt#f- = ■)• Let us define the set of 

generalized curves 

Ce = {(7, p) ■■ 7 € ^’"([0,1], n e ), P £ V, n l 1 ([o, 1])"}. 


Let us consider the following measure Q E on C e 


Q s 



We then have et#Q E = X^fl = g E (t,-) for t £ [0,1]. We define u Q 
and inn® as in ( |14| ) and B respectively, by using test-functions defined 
on Q e . We then have a® = cr E . Indeed, for tp £ C(f2 e ,R d ), we have 


r^=r r 

JJn E Jo 

-j'l 

J 0 J n e 

= <p da E 
in. 


‘ fl(Xt(x )) ■ v E (t,Xt(x))fl{x) dtdx 


p(x) ■ v E (t, x)g E {t, x) dx dt 


which gives the equality. We used the definition of Q E , the fact that 
Xtjj.fl. = g e (t, ■) and that v E g E = a E and Fubini’s theorem. In the same 
way, we have rrfi £ V a . To prove it, we take the same arguments as in 
the end of Step 1 and in the previous calculation. For 1 p £ C(Q. e , R d ), we 
have 



■ v mP (dx, dv) 

N 


g E (t.,X E (x)) 
ip(x) ■ a E (x)dx ) dt 


fl(x)dx J dt 


- £ (/„. 

-['([ „ W <*)) 

Jo \Jn E 

-£(£ 

L 


= ip da E . 


dt 


Moreover, more precisely, we have m jj? (dx,dv) = 5 v .r x \P^(x)dx. Then 
we conclude as in [6]. First for any Lipschitz curve p>, let us denote by ip its 
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constant speed reparameterization, that is, for t G [0,1], (p(t ) = tp(s 1 (t)), 
where 


s(t) = y ——y [ \p(u)\du with Up) = [ \tp(u)\du. 

HW Jo Jo 

For (<p, p) G C, let p be the reparameterization of p i.e. 

Pk(t) '■= ■ | 6 ./g^| t ^i Pfc( s ~ 1 (*)),V* 6 [0,1], k = 1 


Let us denote by Q the push forward of Q through the map {up, p) i-» (ip, p). 
We have and <r® = cr®. Then arguing as in JT5], the L q 

bound on rrfi yields the tightness of the family of Borel measures Q e on 
C([0, l],R d ) x L 1 ([0, 1])^. So Q e *-weakly converges to some measure Q 
(up to a subsequence). Let us remark that Q £ has its total mass equal to 
that of /+, that is, l+e|fi e |. Thus one can show that Q(C) = 1) (due to the 
fact that Q(C) = lim e _ >0+ Q(C e ) = 1). Moreover, we have Q G Q(f~, /+) 
thanks to the *-weak convergence of Q e to Q. Recalling the fact that 
Pjf = rrfi strongly converges in L q to Qk (g G V a ) and due to 

the same semicontinuity argument as in ihi nu, we have m Q (-,Vk (•)) < Qk 
in the sense of measures. Then rrfi (■, u*,(-)) G L q so that Q G Q q (f-, /+). 
It follows from the monotonicity of G(x,v,-) that : 


/ G(x, v, rrfi(x, v)) 9(dx, dv) < / G(x,v, g(x,v)) 9(dx,dv) 

JnxS 11 - 1 J nxS ^- 1 


< inf (121 + 5. 


Letting 5 —> 0 + , we have the desired result. P 

In fact, we showed in the previous proof a stronger result. We proved 
the following equivalence 


Q solves (101 


solves (111 


and moreover, 


(rrfi (• ,Vk(-)))k=i,...,N G V° 


is optimal for (101. We also built a minimizing sequence for (101 from a 
regularization of a solution a of (111 by using Moser’s flow argument. 


3 Characterization of minimizers via anisotropic 
elliptic PDEs 


Here, we study the primal problem (12 I and its dual problem (131. Re¬ 
calling that / = /+—/_ has zero mean, we can reduce the problem (131 
only to zero-mean W 1,V (Q.) functions. Since for (x,v) G O x S d_1 and 
k = G(x,v,-) has a positive derivative at zero, G is strictly 

convex in its last variable then so is Q(x, •) for x G H. Thus Q* is C 1 . 
However Q is not differentiable so that Q*(x,-) is degenerate. By stan¬ 
dard convex duality (Fenchel-Rockafellar’s theorem, see |T5] for instance), 
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we have that min (|12j) = max (131 and we can characterize the optimal 
solution a of (12 I (unique, by strict convexity) as follows 

cr(x) = VG*(x,Vu(x)), 


where u is a solution of (131. In other terms, u is a weak solution of the 
Euler-Lagrange equation 

— div (VC/*(x, Vu(i))) = / in fi, 
V5*(x,Vu(x))-i/f! = 0 on dfl, 

in the sense that 

[ VQ*(x, V«(x)) ■ V<p(x) dx = [ <p(x) df(x), V<p £ W 1,p (fl). 


Let us remark that if u is not unique, a is. 

A typical example is g(x, Vk(x), m) = g k {x, m) = a k (x)m q ~ 1 + 5k with 
5k > 0 and the weights a k are regular and positive. We can explicitly 
compute Q*{x, z). Let us notice that for every x £ fl, z £ R d , we have : 

Q*(x,z)= sup (z ■ u — Q(x, a)) = sup (z ■ cr — inf G(x,q )) 

CTg R d o-6K d e&'PZ 


= sup(« • a - G{x, e)) = sup <V(^' v k (x))Qk - G{x, g) > . 
°,e Q e 1 1 


A direct calculus then gives 


G*(x, z) = V' ( z ■ Vk(x) - 5 k Ck(x)) 

fc! p 


V 

+ > 


where bk = ( akCk) 9-1 . The PDE then becomes 


-EE di [bk(x)vki(x)(Vu ■ Vk{x) - 5kCk(x)) 


p-ii 


= /, 


(15) 


where v k {x) = (v k i(x ),..., v k d.{x)). 

For k = 1,. .., N, Ql(x,z) = bk ^ {z-Vk(x) — 5k) p + vanishes if z-v k {x) £ 
] — oo, 5kC k {x)] so that any u whose the gradient satisfies Vu(x) ■ v k {x) £ 
] — oo, 5kCk(x)], Vx £ fi, k = 1, ..., N is a solution of the previous PDE 
with / = 0. In consequence, we cannot hope to obtain estimates on 


the second derivatives of u or even oscillation estimates on Vu from 1151. 


Nevertheless we will see that we have some regularity results on the vector 
field a = ( 01 ,... ,aa) that solves ( |12| ) in the case where the directions and 
the volume coefficients are constant, that is, 


o-(x) = ^2 [6fc(x)(V«(x) • v k - 5 k Ck) p + x ] v k , 


for every x £ SI. 
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4 Regularity when the v^s and Ck s are 
constant 


Our aim here is to get some regularity results in the case where the v k s 
and the c k s are constant. We will strongly base on [6] to prove this 
regularity result. Let us consider the model equation 

N 

- ^ div ((Vb(i) ■ v k - SkCkY+^Vk) = f, (16) 

fc = l 

where Vk € § d_1 , Ck > 0 and bk = 1 for k = 1,..., N. Define for z £ 


N 

F{z) = ^ Fk(z), with Fk(z) = (z • v k - £fcCfc)+ _1 ^fc (17) 

k=1 

and 

N 

H{z) = y ^H k (z), with H k (z) = (z ■ v k - S k c k )lv k - (18) 

fc=1 

Here we assume only p > 2. We have the following lemma that establishes 
some connections between F and H. 

Lemma 4.1. Let F and G be defined as above with p > 2, then for every 
(z, w)€R d x R d , the following inequalities are true for k = 1,..., N 

\F k (z)\ < |*r\ (19) 

\F k (z) - F k (w)\ < (p-1) (\H k {z)\ V -^ + \H k {z)\^ | H k (z) - H k (w)\, 

( 20 ) 

and 

(. F k (z) - F k (w)) -(z-w)> 4l H k (z) - H k (w )| 2 . (21) 

pZ 

Proof. The first one is trivial. For the second one, from m one has the 
general result: for all (a, 6) £ R d x R d , the following inequality holds 


|a| p_2 a — |6| p_2 fe| < (p-1) (\a\^ + |6|V) ||a|^o - |6|^&| . (22) 


Choosing a = (z • v k — 5 k c k )+v k and b = (w • v k — 5 k c k )+v k in (J22J), we 
then obtain (201. 

Let us now prove the third inequality. It is trivial if both z ■ v k and 
w ■ v k are less than 5 k c k . If z ■ v k > 5 k c k and w ■ v k < 5 k c k , we have 

(■ F k (z)-F k (w))-(z-w ) = (z-Vk-SkCk) 1 ^ 1 (z-v k -w-v k ) > ( z-v k -6 k c k )+ = \H k (z)\ 2 . 


For the case z-v k > 8 k c k and w-v k > 5 k c k , we use the following inequality 
(again [17] ') 

(|a| p_2 a - |b| p_2 &) • (a - 6) > 4 ^al^a - \b\^b^ . 
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Again taking a = (z ■ v k — 8kC k )+v k and b = (w ■ Vk — 8kCk)+Vk, we have 
that 

\\H k (z) - H k {w)\ 2 

pZ 

< (\F k (z)\ - |F fe (w)|)u fe • {(z • Vk - 5 k Ck)+ - (ui ■ v k - S k c k )+)v k 
= (\F k (z)\ - |F fc (w)|)(z - w) ■ v k , 

which gives ( |2 1| ) . □ 

Let us fix / € where q is the conjugate exponent of p and let 

us consider the equation 

— div_F(Vu) = /. (23) 


Thanks to Nirenberg’s method of incremental ratios, we then have the 
following result that is strongly inspired of Theorem 4.1 in [6]: 

Theorem 4.1. Let u £ W, l m p (12) be a local weak solution of (|23|). Then 


H := L/(Vu) £ (f l). More precisely, for every k = 1, ... ,N,TL k '■= 

Hk (Vit) £ <’ 2 (fi). 

Proof. For the sake of clarity, write T := F’(Vu) and similarly, Fk,Hk 
(note that T k £ L q oc (Ll) and Hk £ L 2 oc (Pl) due to (19 H 201. Let us define 
the translate of the function i p by the vector h by th<P := y?(- + h). Let 
ip £ W 1 ' q { fl) be compactly supported in Q. and h £ R d \{0} be such that 
\h\ < dist(supp(i^),R d \{0}), we then have 


r h T - T 
\h\ 


Vpdx = 


/ 

Jo. 


Thf ~ f 


pdx. 


(24) 


Let u> <s ojo <s fl and £ £ C£°(C2) such that supp(£) C wo,0 < £ < 1 
and £ = 1 on uJ and h £ R d \{0} such that \h\ < ro < |dist(u>o,R d \f2). 
In what follows, we denote by C a nonnegative constant that does not 
depend on h but may change from one line to another. We then introduce 
the test function 

<p = £ 2 \h\~ 1 (r h u - u), 

in (241. Let us fix ui' := u>o + B(0,ro). It follows from u £ W^(f(f2),/ £ 
Wtm an d the Holder inequality that 


\h\~ 


'■ f (T h F-F)-(e(r h \/u-\/u) + 2^ar h u-u)) < ||V/||z,, (w0 1 |V«|Up( u 
Jn 


The left-hand side of the previous inequality is the sum of 2N terms 
In + I 12 + ■ ■ ■ + Ini + In 2 where for every k = 1 ,..., N, 


hi ■= H 


7 

Jn 


f{F k {T h Vu) - F k (Vu) ■ (rhVu - Vu), 


and 


Ik 2 ■= \h\~ 2 [ f(F k (T h Vu) - F k (Vu) ■ Vftfou - u). 

Jn 

Let k = 1,..., N fixed. We will find estimations on hi and I k 2 - Due to 
(201, hi satisfies: 

hi > ^mh\~\r h H k -Hk)\\l2. 

pZ 
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For Ik 2 , if p > 2, it follows from (211 and the Holder inequality with 
exponents 2 ,p and 2 p/ ( p — 2) that 


|/*a| < \ h \~ 2 ^ \^^\\ThU-u\\T h Hk-Hk\ {\r h H k \^ +\H k \^ 
<C\\\h\- 1 {T h u-u)\\ LP ^ o) \\m~ 1 (.T h Uk-U k )\\ L i ( [ \H k \ 2 + \r h H k \ 

\J cuo 

<GU\h\~\r h H k -Hk)\\ L ^ 

and if p = 2, we simply use Cauchy-Schwarz inequality and we get : 
\h2\<CU\h\- 1 (TH'H k -'H k )\\ L 2. 


p-2 

2 P 


Bringing together all estimates, we then obtain 


E 


ThHk — Hk 


k =1 

and we finally get 


h 




< c i+E 


Th'Hk — 'Hk 
h 



E 


Th'Hk — Hk 
h 


2 


<C, 


for some constant C that depends onp, ||/||vyi>«> and the distance 

between cj and <9S2, but not on h. We have the desired result, that is, 
Hk £ W'jo’c (H), for k = 1,..., N, and so H also. □ 

If we consider the variational problem of Beckmann type 

inf If inf c k (-g q k + S k Qk] ■ -div a = f 1 , (25) 

^6i=(Q) \q ) 


we then have the following Sobolev regularity result for the unique mini- 
mizer that generalizes Corollary 4.3 in [6]. 

Corollary 4.1. The solution a of ( |25| ) is in the Sobolev space 
where 


' 2 

any value < 2, 
dp 

. dp — (d + p) + 2 ’ 


ifp = 2 , 

if p > 2 and d = 2, 
if p > 2 and d > 2. 


Proof. By duality, we know the relation between a and any solution of 
the dual problem u 


N 

a = ^(Vu • v k - 5 k Ck) p i r 1 v k ■ 

k =1 

Since u £ W 1 ’ q (Pl) is a weak solution of the Euler-Lagrange equation (|16|), 
using [Theorem 4~T[ and |Lemma i~T| we have that the vector fields 

P 

H k {x) = (Vu(x) • v k - S k Ck)lv k , k = 1,... ,N, 
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are in W’ 1 *’ < ?(12). We then notice that a = ak with 

oh = \Hk\^Hk, k = l,...,N. 

The first case is trivial: we simply have Ok = Hk £ For the other 

cases, we use the Sobolev theorem. If p > 2 and d > 2 then Hk £ L\* c m 
with 

1 _ 1 1 
2* ~ 2 ~ N' 


Applying (201 with 2 = r^Vu and w = Vu, we have 

< (v — 


'l~h&k &k 


ip- 1) (\T h 'H k \ F r +| Hk\” p ^ 


Th'Hk — Hk 


\h\ 


2*P 

- p — 2 i 


Since [Hk\ p £ L p oc (fl), we have that the right-hand side term is in 
Lr oc (n) with r given by 

1 _ p-2 1 

r ~ 2 *p + 2' 

We can then control this integral 


ThOk &k 


\h\ 


dx. 


For the case p > 2 and d = 2, it follows from the same theorem that 
Hk £ Tf oc (12) for every s < +00 and the same reasoning allows us to 
conclude. □ 

This Sobolev regularity result can be extended to equations with weights 
such as 


- ^2 div (bk(x){’Vu(x) ■ vk - 5kCk) p + 1 Vk) = /■ 


(26) 


An open problem is to investigate if one can generalize this Sobolev 
regularity result to the case where the Vk s and Ck’s are in C ,1 (12). 


5 Numerical simulations 

5.1 Description of the algorithm 

We numerically approximate by finite elements solutions of the following 
minimization problem: 


inf J(u):=G*(Vu) -(/,«> (27) 

uEW^-’P (f2) 

with G*(4>) = f n G*(x, $(*)) dx for $ 6 L v (Q.) d and ( f,w ) = f n u df 
for w £ L p (12). Let us recall that 12 is a bounded domain of R d with 
Lipschitz boundary and / = /+ — /_ is in the dual of W 1,p ( 12) with zero 
mean f Q f = 0. We will use the augmented Lagrangian method described 
in |5] (that we will recall later). ALG2 is a particular case of the Douglas- 
Rachford splitting method for the sum of two nonlinear operators (see [T8j 
or more recently [20]). ALG2 was used for transport problems for the first 
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time in [4], Let a regular triangulation of f 1 with typical meshsize h, let 
Eh C W 1,p (£l) be the corresponding finite-dimensional space of P 2 finite 
elements of order 2 whose generic elements are denoted Uh- Moreover, we 
approximate the terms / by fh E Eh (again with ( fh , 1} = 0) and G by a 
convex function G^. Let us consider the approximating problem 


inf Jh{u h ) '■= Gl(Vuh) - {fh,Uh)- (28) 

v-hSE h 


and its dual 


sup {-G h {(Th) ■ —di v h (a h ) = fh} (29) 

where Fh is the space of Pi finite elements of order 1 and —div/ l (cr/ l ) may 
be understood as 


{a h .yu h ) F d = -(di v h {a h ),u h )E h - 


Theorem 5.1. If Uh solves ( |28| ) then up to a subsequence, Uh converges 
weakly in W 1,p ( fl) to a u that solves (|27[) as h —> 0. 


It is a direct application of a general theorem (see [5] and [14t for similar 
results and more details). Using the discretization by finite elements, (271 
becomes 


inf J(u) :=F(u) + G*(A u) 

u(zM. n 


(30) 


where F : R n —> R U {+ 00 }, G : R m —> R U {+ 00 } are two convex l.s.c. 
and proper functions and A is an m x n matrix with real entries. A is the 
discrete analogue of V. The dual of (|30[) then reads as 


sup — F*(—A T cr) — G(cr) (31) 


' x R m satisfies the primal-dual extremal- 


We say that a pair (u, cr) 
ity relations if: 

— \ T a £ <9F(u), a G 9G*(Au). (32) 

It means that u solves (301 and that 0 solves (311 and moreover, (|30| 


and (311 have the same value (no duality gap). It is equivalent to find a 
saddle-point of the augmented Lagrangian function for r > 0 (see 
for example) 


L r (u,q,a) := F(u)+G*(g)-|-CT-(A'u— q) + l~\Au—q\ 2 ,V(u,q,a) eR"xR”xR ra . 

(33) 

It is the discrete formulation of the corresponding augmented Lagrangian 
function 


L r (u,q,a) := / g*(x, q{x)) dx - (u, f) + (a, Vm - q) + T -\Vu(x) - q(x)\ 2 

J f2 A 

(34) 

and the variational problem of ( |30| ) is 

inf < [ G*(x, q(x)) dx — [ u{x)f{x) dx 1 . (35) 

u ’i U n Jn J 

subject to the constraint that Vu = q. 

The augmented Lagrangian algorithm ALG2 involves building a se¬ 
quence (u k , q k ,a k ) G R x R d x R d from initial data ( u°,q° , cr°) as follows: 
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1. Minimization problem with respect to u: 

u k+1 p= argmin ueH „ |f(m) + a k ■ An + ^|Vm - <2*j 2 )} 

That is equivalent to solve the variational formulation of Laplace 
equation 

—r(Vu fe+1 — div(g fc )) = / + div(<r fe ) in D. 
with the Neumann boundary condition 

du k+1 


dv 


k k on 

= rq • v — a • u on oil. 


This is where we use the Galerkin discretization by finite elements. 

2. Minimization problem with respect to q\ 

q k+1 := argmin 9SRd |G*(g) - o k ■q+ ^|VM fc+1 - q\ 2 ) j 

3. Using the gradient ascent formula for a 

fc + 1 k , fe + 1 fc + l\ 

a =a + r(VM — q ). 

Theorem 5.2. Given r > 0. If there exists a solution to the primal- 


dual extremality relations (321 and A has full column-rank then there 
exists an (u,d) G R n x R m satisfying (321 such that the sequence 


( u k ,q ,u k ) generated by the ALG2-scheme above satisfies 

u k —> u,q k —> Am, a k —¥ d as k —^ +oo. (36) 

We directly apply a general theorem whose proof can be found in m 
(Theorem 8), following contributions of [13, 14, 48] to the analysis of 
splitting methods. 

5.2 Numerical schemes and convergence study 

We use the software FreeFem+T (see [T6]) to implement the numerical 
scheme. We take the Lagrangian finite elements and notations used in 


Subsection 5.1 P 2 FE for mj, and Pi FE for ( qh, o^)- A.Uh is the projection 


on Pi of the operator A, that is, Vuh- The first step and the third one 
are always the same and only the second one varies with our different 
test cases. We indicate the numerical convergence of ALG2 iterations by 
the • k superscript and the convergence of finite elements discretization by 
the -h subscript. For our numerical simulations, we work with the space 
dimension d = 2 and we choose for Q a 2D square (x = (* 1 , 2 : 2 ) G [0, l] 2 ). 
We make tests with different / : 

f 1 ._ e -M*«x 1 -0.75) 2 +(x2~0.25) 2 and j-1 e -40*((®i -0.25) 2 + (x 2 -0.65) 2 ) 


j2 e -40*(( a ;i-0.5) 2 + ( a;2 -0.15) 2 ) and , = g —40«((a: 1 -0.5) 2 + (x 2 -0.75) 2 ) 
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In the third case, we take /£ a constant density and f 'l is the sum of 
three concentrated Gaussians 

f 3 {x 1,X 2 ) = 400 *d a:— 0-25) 2 + (y—0.75) 2 ) + e -400»((ic-0.35) 2 +(!/-0.15) 2 ) 

+ e -400«((:c-0.85) 2 + (i/-0.7) 2 ) 

We also make tests with non-constant Ck : 

g\x 1 ,x 2 ) = 3 - 2 * e -io*((x 1 -o.5) 2 +(!/2 -o.5) 2 )_ 

As specified above, we use a triangulation of the unit square with 
n = 1/h element on each side. We use the following convergence criteria: 

1. DIV.Error = (/„ h (diver* + /) 2 ) ^ is the L 2 error on the divergence 
constraint. 

2. BND.Error = (/ enfe (er£ • i/) 2 ) 1/2 is the L 2 (dCi h ) error on the Neu¬ 
mann boundary condition. 

3. DUAL.Error = max^ \G(xj, + G*{xj, V«h(%')) — V u^ixj ) ■ 

a h( x j)\ where the maximum is with respect to the vertices Xj. 

The first two criteria represent the optimality conditions for the mini¬ 
mization of the Lagrangian with respect to u and the third one is for 
maximization with respect to cr. 

We make tests for two models. In the first one, the directions are 
the same as in the cartesian model and the volume coefficients are not 
necessarily constant. In the second one, the directions are the same than 
in the hexagonal one and the volume coefficients are equal to 1 (it is 
simpler to compute Q(x,a)). That is, Vk = exp(ifc-7r/3) and SkCk = 1 for 
k = 1,... ,6. We call these models still the cartesian one, the hexagonal 
one respectively. The cartesian one is much easier since we can separate 
variables. G = Gi + G 2 with G i(x,q) = ^{\qi\ — <5iCi(ar))fJ_ so that the 
second step of ALG2 is equivalent to solve the pointwise problem 

inf -(k|-c(*))+ + ^l'7-9 fc i 2 

q p Z 

where q k = Vu k+1 + This amounts to set q k+1 = A q k and to solve 
this equation in A 

(X\q k \-c(x)) p + - 1 +rX\q k \=r\q k \=0 

with A > 0. We can use the dichotomy algorithm. 

For the hexagonal one, we use Newton’s method. Since the function 
of which we seek the minimizer has its Hessian matrix that is definite 
positive, we can use the inverse of this Hessian matrix. 

We show the results of numerical simulations after 200 iterations for 
both models. 

We notice that length of arrows are proportional to transport density. 
Level curves correspond to the density term of the source/sink data to 
be transported. In|Figure 3[ the case p = 1.01 means that there is much 
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Test case 

DIV.Error 

BND.Error 

DUAL.Error 

Time execution (seconds) 

1 

8.4745e-05 

0 

3.6126e-06 

436 

2 

2.2536e-05 

8.8705e-04 

3.0663e-05 

4764 

3 

5.2141e-05 

1.4736e-04 

1.1556e-02 

792 

4 

1.1823e-05 

7.6776e-04 

8.7412e-06 

170 

5 

1.1629e-05 

0 

9.7498e-04 

285 

6 

3.1544e-04 

1.0958 

7.8350e-07 

445 

7 

4.1373e-04 

1.1710 

4.8113e-04 

4657 


Table 1: Convergence of the finite element discretization for all test cases. 


congestion. The case p = 2 is reasonable congestion and in the last one 
p = 100, there is little congestion. When there are obstacles, the criteria 
BND.Error is not very good. Indeed, the flow comes right on the obstacle 
and it turns fast. In the other side of the obstacle, the flow is tangent to 
the border. Many other cases may of course be examined (other boundary 
conditions, obstacles, coefficients depending on x, different exponents p for 
the different components of the flow...). 

Acknowledgements The author would like to thank Guillaume Car- 
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_V«c Vilu< 

■ 

■3 0620872 


■3.248349 
■3.310436 
■0 372523 


§3.682959 
13.745046 
13.807133 
(1.86922 
13.931307 
13.993395 
jl .05548 
|1. 11757 
|1. 17966 



_ Vec Vilu* 

■' 

■3 124866 

■> 249-,-32 


■3.499463 

■3.624329 

■'749195 


■137352 

■1.49839 
■1.62326 
■1.74812 
■1.87299 
■1.99785 
■2.12272 
■2 24758 
■2.37245 



Figure 3: Test cases 3, 4 and 5: cartesian case (d = 2) with / = / 2 , Ck constant 
and p = 1.01, 2,100. 
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^Vec Value 

■1 823711 


■2 47113 
■3.29484 
■4.11856 
■4.94227 

■f. 765.9ft 


■6.23711 

■ .06082 

■>.88453 

■10.7082 

■11.532 

■12.3557 

■13.1794 

■14.0031 

■14.8268 

■15.6505 


Figure 4: Test case 6 : cartesian case (d = 2) with / = f 1 , c*, = g 2 , p = 3 and 
two obstacles. 



^Vec Value 

■l 06349 


■ 19046 

■ 25395 

■ .31744 

■ 38093 


■11.6984 

■12.7619 

■13.8253 

■14.8888 

■15.9523 

■17.0158 

■18.0793 

■19.1428 

■0.2063 


Figure 5: Test case 7 : hexagonal case (d = 2) with / = f 1 , Ck constant, p = 3 
and an obstacle. 
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